program lusolve_lapack
	implicit none
	integer, parameter :: N = 500
	integer :: i, j, k
	integer,allocatable :: ipiv(:)
	integer :: info
	double precision, allocatable :: A(:, :), b(:), x(:), Abak(:, :)
	integer :: scc1, scc2, sccm, sccr 
	double precision :: dtemp
	real :: rtemp

	call system_clock( scc1, sccm, sccr )

	! Initialize
	allocate( A(N, N) )
	allocate( Abak(N, N) )
	allocate( b(N) )
	allocate( x(N) )
	allocate( ipiv(N) )
	
	do i = 1, N
		do j = 1, N
			call random_number(rtemp)
			A(i,j) = dble(rtemp)
		end do
		A(i,i) = A(i,i) + 10d0
		call random_number(rtemp)
		b(i) = dble(rtemp)
	end do

	Abak = A
	! Solve
	call DCOPY( N, b, 1, x, 1 )
	call dgesv( N, 1, A, N, ipiv, b, N, info )
	call DSWAP( N, b, 1, x, 1 )
	call DGEMV( 'n', N, N, 1d0, Abak, N, x, 1, -1d0, b, 1 )

	write(*,*) "Error:"
	write(*,*) maxval( dabs(b) )
	call system_clock( scc2, sccm, sccr )
	write(*,*) "Time:"
	write(*,*) dble(scc2 - scc1) / dble(sccm)

	deallocate( A )
	deallocate( b )
	deallocate( Abak )
	deallocate( x )
	deallocate( ipiv )

end program

